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■ ABSTRACT 

: 

O ; 

^ . We use numerical simulations to study the effect of nonlinear MHD waves 

\0 . in a stratified, self-gravitating molecular cloud that is bounded by a hot and 

^ , tenuous external medium. In a previous paper, we had shown the details of a 

p . standard model and studied the effect of varying the dimensionless amplitude 

of sinusoidal driving. In this paper, we present the results of varying two other 
important free parameters: /3o, the initial ratio of gas to magnetic pressure at the 
! cloud midplane, and z/q, the dimensionless frequency of driving. Furthermore, we 

J> ! present the case of a temporally random driving force. Our results demonstrate 

^ , that a very important consideration for the actual level of turbulent support 

against gravity is the ratio of driving wavelength Aq to the the size of the initial 
non-turbulent cloud; maximum cloud expansion is achieved when this ratio is 
close to unity. All of our models yield the following basic results: (1) the cloud 
is lifted up by the pressure of nonlinear MHD waves and reaches a steady-state 
characterized by oscillations about a new time-averaged equilibrium state; (2) 
after turbulent driving is discontinued, the turbulent energy dissipates within a 
few sound crossing times of the expanded cloud; (3) the line-width-size relation 
is obtained by an ensemble of clouds with different free parameters and thereby 
differing time-averaged self-gravitational equilibrium states. The best consistency 
with the observational correlation of magnetic field strength, turbulent line width. 
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and density is achieved by cloud models with /3o ~ 1- We also calculate the 
spatial power spectra of the turbulent clouds, and show that significant power is 
developed on scales larger than the scale length Hq of the initial cloud, even if the 
input wavelength of turbulence Aq ~ -f^O) The cloud stratification and resulting 
increase of Alfven speed toward the cloud edge allows for a transfer of energy 
to wavelengths significantly larger than Aq. This explains why the relevant time 
scale for turbulent dissipation is the crossing time over the cloud scale rather 
than the crossing time over the driving scale. 

Subject headings: ISM: clouds — ISM: magnetic fields — methods: numerical — 
MHD — turbulence — waves 



1. Introduction 

Interstellar molecular clouds have nonthermal line widths which are typically supersonic 
and increase with the scale of the cloud size that is measured (e.g., Larson 1981; Myers 
1983; Solomon et al. 1987; Brunt & Heyer 2002). For a cloud scale ~ 1 pc, the velocity 
dispersion a 1 km s~^, which is several times greater than the isothermal sound speed Cg = 
0.19 km s^^ for a typical cloud temperature T = 10 K (Goldsmith & Langer 1978). When 
detectable through the Zeeman effect, magnetic field strengths arc such that their infiucnce 
is comparable to that of gravity (Crutchcr 1999) and velocity dispersions are correlated with 
the mean Alfven speeds Va such that a /Va ~ 0.5 (Crutcher 1999; Basu 2000). Furthermore, 
maps of polarized emission from dust imply that the magnetic field morphology is relatively 
well-ordered and therefore not dominated by turbulence (e.g., Schleuning 1998; Matthews 
& Wilson 2000; Houde et al. 2004). It seems possible that turbulence within molecular 
clouds comprises a spectrum of nonlinear magnctohydrodynamic (MHD) waves generated 
from various localized sources and from self-gravitational motions. In contrast to the large- 
scale structuring by turbulence and magnetic fields, we note that the embedded dense cores 
in regions of low-mass star formation have characteristically subsonic or transonic turbulence; 
evidence for gravitational collapse in cores also comes in the form of subsonic infall motions 
where detectable (see discussion and references in Myers 2005). 

Although not possessing internal motions as highly supersonic as some larger, more 
distant molecular clouds, the well-resolved Taurus molecular cloud complex (distance ^ 140 
pc) reveals a hierarchy of density and size scales that is characteristic of star-forming regions 
(see Ungerechts & Thaddeus 1987; Mizuno et al. 1995; Onishi et al. 1998; Onishi et 
al. 2002). The cloud has a common (turbulent) envelope of molecular column density 
N '^Ix 10^^ cm~^ — 2 X 10^^ cm^^ (corresponding to visual extinction Ay ~ 1 — 2), embedded 
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parsec-scale dark clouds with Ay > 3 and velocity dispersion a 0.6 km s~^ ^ 3Cs, which 
in turn contain dense regions with an average Ay ~ 5 that are forming small clusters of 
stars. The typical separation of young stellar objects within the dense regions may be 
understood by a simple Jeans-like fragmentation process for an isothermal dense sheet or 
filament (Hartmann 2002). Such an approach ignores the effect of turbulence, which is not 
measured to be large in those regions anyway. However, the subsonic infall motions imply 
that a fragmentation process mediated by magnetic fields and ion-neutral friction may be 
more appropriate (Basu & Ciolek 2004). 

In this paper, we extend a previous set of models (see Kudoh & Basu 2003, hereafter Pa- 
per I) which approximate conditions within the molecular cloud envelopes and outer regions 
of dark clouds. These regions are expected to be sufficiently ionized by the far-ultraviolet 
background starlight that the ideal MHD limit of fiux-freezing is applicable (McKee 1989). 
Some key questions about cloud envelopes are as follows. What is the rate of dissipation 
of the waves and what is the primary dissipation mechanism? Is continual driving of the 
turbulence necessary to explain their ubiquitous presence in all but the densest regions? 
Furthermore, observed properties of the turbulence such as the line- width-size relation (e.g., 
Solomon et al. 1987) and the correlation between magnetic field strength, line width, and 
density (Myers & Goodman 1988; Basu 2000) deserve an explanation. 

Since the late 1990's, a set of three-dimensional models of MHD turbulence have been 
developed which use a local approach, i.e., use periodic boundary conditions (e.g.. Stone, 
Gammie, & Ostriker 1998; Mac Low ct al. 1998; Mac Low 1999; Padoan & Nordlund 1999; 
Ostriker, Stone, & Gammie 2001). This means that the modeled regions represent a very 
small region deep within a molecular cloud. However, they are driven with highly supersonic 
motions characteristic of large-scale molecular cloud structure. Although these models are 
local, the turbulent driving is global, in that large amplitude momentum fiuctuations are 
input at every grid point. The fiuctuations are drawn from a Gaussian random field with 
a specified root mean squared amplitude and power spectrum. Although this method is 
convenient, a more likely scenario for molecular cloud turbulence is that it is driven from 
local sources and propagates to fill the clouds. A main result of the periodic-box models 
is that the decay time of the MHD turbulence is comparable to the crossing time over the 
driving scale of the turbulence. For instance. Stone et al. (1998) find a turbulent dissipation 
rate Fturb ~ PoO""^/Ao, where po is the fixed mean density in the periodic simulation box, a is 
the one-dimensional velocity dispersion, and Aq is the driving scale of turbulence. A direct 
application of this result to large molecular clouds has been questioned by Basu & Murah 
(2001), who point out that dissipation proportional to the crossing time across internal 
driving scales would lead to a CO emission luminosity Leo far exceeding that observed for 
clouds whose size is much larger than the driving scale. Alternatively, they point out that 
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the observed scaling of Leo with cloud mass can be explained if the relevant length scale for 
dissipation is the largest scale L of any cloud; this also yields the correct order of magnitude 
for Leo itself. In the context of a periodic box model, this means that the driving would have 
to occur only on the largest scales of each cloud. McKcc (1999) also questions extremely fast 
dissipation, on the grounds that: (1) turbulence should not be so ubiquitous (and present 
for example in clouds with or without embedded OB associations, a possible major source of 
stellar-driven turbulence) if it is decaying so rapidly; (2) dissipation times of at least a few 
cloud crossing times (as opposed to a free-fall time) may explain the relatively low Galactic 
star formation rate; and (3) current simulations cannot resolve the turbulence within the 
dense clumps that quickly develop. Sugimoto, Hanawa, & Fukuda (2004) have supported 
the last point with second-order accurate three-dimensional numerical studies of Alfven 
wave propagation, finding that ~ 30 grid points per wavelength are necessary to reduce 
dissipation to 1% per wavelength propagated. Such a high resolution remains prohibitive 
for three-dimensional simulations that seek to model large-scale cloud dynamics and dense 
core formation. Cho & Lazarian (2003) have also questioned a common interpretation of 
three-dimensional periodic box models, i.e., that the fast dissipation is due to the coupling 
of Alfven modes to the more dissipative fast MHD and slow MHD modes (e.g.. Stone et al. 
1998). Through modal analysis of their three-dimensional simulations, they find that the 
coupling of modes is rather weak in the inertial range of turbulence and that dissipation of 
Alfven waves is dominated by a nonlinear cascade. Furthermore, Cho, Lazarian, & Vishniac 
(2002) have found that in the case of incompressible MHD turbulence, a directionality in the 
flux of energy (as may be expected with localized driving sources) leads to a significantly 
lower dissipation rate. 

In addition to the detailed analysis of individual wave modes as described above, new 
approaches to the numerical study of MHD turbulence include extremely high resolution 
models of the global structure of molecular clouds (Kudoh & Basu 2003; Folini, Heyvaerts, 
& Walder 2004). These one-dimensional models employ a numerical resolution that remains 
beyond the reach of three-dimensional simulations. Kudoh & Basu (2003, hereafter Paper 
I) have modeled the global structure of a one-dimensional cold cloud that is embedded in 
a hotter external medium. The main advantage over a local (periodic box) approach is the 
inclusion of the effects of cloud stratification and a cloud edge. Furthermore, the turbulent 
driving is localized and the propagation of disturbances within the cloud as well as their 
leakage from the cloud can be studied. Our model in Paper I had 50 numerical cells within 
the scale length Hq of the initial state, and 70 cells within the wavelength Aq associated with 
turbulent driving. Results included the generation of significant density structure and shock 
fronts within the cloud; however, most of the wave energy remained in transverse, rather than 
longitudinal motions (in qualitative agreement with the results of Cho & Lazarian 2003) and 
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the time-averaged density profile of the cloud (which may be similar to an observed spatial 
average of many layers along the line-of-sight) was relatively smooth. The localized driving 
led to turbulence which quickly filled the cloud, causing cloud expansion and large-scale 
oscillations in the outer regions; all this is aided by a clear transfer of energy to the largest 
scales of the expanded cloud, due to the stratification and consequent increase of Alfven 
speed toward the cloud edge. The work of Folini et al. (2004), who model a very similar 
system with similar turbulent driving, confirm the existence of highly disturbed density 
structure, oscillatory motions, and a dominance of wave energy in transverse modes. They 
focus most extensively on the density structuring and demonstrate that progressively finer 
scale structuring and enhanced cloud support is revealed by higher resolution simulations; 
for typical dimensional values of physical parameters, their numerical resolution is 0.001 pc, 
as in our models in Paper I. In Paper I, we also considered many global properties of the 
clouds. We found that random motions were greatest in the outer, low-density regions of 
the clouds in a way that the velocity dispersion a and cloud scale L satisfied the observed 
line-width-size relation. Our model also showed that even when turbulence decays away due 
to lack of driving, the large-scale oscillations are the longest lived component. These results 
can help to explain the observations of apparent oscillations of B68 (Lada et al. 2003), a dark 
cloud in which nonthermal motions are small. We also showed that the overall dissipation 
of turbulent energy (after driving is discontinued) is proportional to the crossing time across 
the expanded cloud, rather than the crossing time for the driving wavelength Aq. This result 
points to a way out of the problems with internal driving pointed out by Basu & Murali 
(2001). Stratification leads to turbulence at the largest scales, so that the cloud scale L 
is more relevant to turbulent dissipation than the internal driving scale Aq, as needed on 
empirical grounds. We quantify the effect of turbulence on large scales further in this paper 
by considering the power spectrum of fluctuations in the turbulent clouds. 

This paper expands upon the study of Paper I and summarizes the results of a parameter 
study. The numerical model is presented in § 2. The results of the simulations are described 
in § 3. A summary and discussion are presented in § 4. 

2. The Numerical Model 

2.1. Overall Approach 

The basic model is the same as presented in Paper I. We consider a molecular cloud 
that is threaded by a large-scale magnetic field and assume ideal MHD motions. We assume 
a driving force near the midplane of the cloud and follow the dynamical evolution of the 
vertical structure of the cloud. Our model can be characterized as 1.5-dimensional, since 
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we calculate variations in only one direction but allow vector quantities to have both a 
longitudinal and a transverse component (see details in § 2.2). We also assume isothermality 
for each Lagrangian fluid element. This means that the temperature does not change in time 
for each fluid element as it moves through Eulerian space. Since we have a two-temperature 
model (cloud and intercloud medium), this is not the same as a purely isothermal model. 



2.2. Basic Equations and Physical Variables 

Wc use local Cartesian coordinates {x, y, z) on the molecular cloud, where we set z 
to be the direction of the large-scale magnetic field. According to the symmetry of the 
1.5-dimensional approximation, we set 

dx dy ^ ^ ^ 

The above symmetry and the divergence-free condition on the magnetic field imply 

— constant, (2) 

where is the z-component of the magnetic field that threads the molecular cloud. More- 
over, from the assumption of linear polarization of the waves, we can set 

v,^B, = Q (3) 

without loss of generality, where and B^ are the x-components of the velocity and magnetic 
field, respectively. 

Under this assumption, we solve the MHD equations numerically. In the equations, t is 
the time, G is the gravitational constant, k is Boltzmann's constant, m is the mean molecular 
mass, p is the density, P is the pressure, T is the temperature, is the 2;-component of the 
gravitational field, is the 2;-component of velocity, Vy is the y-component of velocity, and 
By is the y-component of the magnetic field. In the equations, the energy equation is assumed 
to be 

which quantifies the assumption of isothermality for each Lagrangian fiuid element. 
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2.3. Initial Conditions 

As an initial condition, we assume hydrostatic equilibrium of a self-gravitating one- 
dimensional cloud. The hydrostatic equilibrium is calculated from the equations 

p dz 



and 



subject to the boundary conditions 



t = (6) 



and 



g,{z = 0) = 0, (8) 
p(z = 0) = po, (9) 
kT 

P{z = 0) = Po — , (10) 

m 

where po and Tq are the initial density and temperature at ^ = 0, respectively. 

In order to solve the above equations, we need to assume an initial temperature distri- 
bution. If the temperature is uniform throughout the region, we have the following analytic 
solution ps found by Spitzer (1942): 

Psiz) = posedi^iz/Ho), (11) 



where 



is the scale length, and 



Ho = , (12) 



CsQ = \ — (13) 

V m 

is the isothermal sound speed for temperature To- 

However, an isothermal molecular cloud is usually surrounded by warm or hot material, 
such as neutral hydrogen or ionized gas. Therefore, we assume the initial temperature 
distribution to be 



r(^) = ro + ^(r,-ro) 



1 -|- tanh 



(14) 
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where we take Tc — 100 To, -^c = 3ifo, and Zd — 0.2Hq throughout the paper. This distribu- 
tion shows that the temperature is uniform and equal to Tq in the region of < ^ < = 3-ffo 
and smoothly increases to another uniform value Tc = IOOTq at z ~ = ^Hq. By using this 
temperature distribution, we can solve the ordinary differential equations (5)- (7) numerically. 

We also assume the following initial conditions: 

v^{z) = Vy{z) = 0, (15) 

By{z) = 0, (16) 

B,{z)=Bo, (17) 

where Bq is a constant. According to equation (2), is spatially uniform and independent 
of time throughout the calculations. 



2.4. Driving Force 

We introduce a perturbation into the initially hydrostatic cloud by adding a driving 
force, F{z,t), into the y-component of the momentum equation as follows: 

In this paper, wc use two forms of the driving force. The first one, which is the same as 
in Paper I, is the sinusoidal driving force defined by 

pad(^) sin(27rz/ot) exp[-(^)2] {t < lOto) 
F{z,t)^{ padSm{2nuot)exp[-{fJ^] " {lOto < t < AOto) (19) 
" (i > 40io), 



and 



to = Ho/cso. (20) 



In this form, is the amplitude of the induced acceleration, uq is the frequency of the driving 
force, and Za represents the region in which we input the driving force. 

The other form we use is the random driving force: 

pa,(^)ran(t) exp[-(^)2] (t < lOto) 
F{z, t) = { paars^nit) exp[-(^)2] (lO^o < t < AOto) (21) 

" (t> AOto), 
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where ran(i) is a random number between -1 and 1 which is determined at each time step. 

These equations show that we input the driving force near the midplane of the cloud, 
and increase the maximum driving force hnearly with time until t — lOtg, and maintain it 
to be constant during lO^o <t< AOto. After t — 40io, we terminate the driving force. 



2.5. Boundairy Conditions and Numerical Technique 

We used a mirror symmetric boundary condition at ^ = and a free boundary at 
z — ^out, the outer boundary of the calculation. In order to remove the reflection of waves 
at the outer boundary, we set Zo^t to be a large value, and use nonuniform grid spacing 
for large z. For z < SOHq, the grid spacing is Azi — 0.02i7o; most of the 4000+ points in 
a typical simulation are concentrated in this region. Such high resolution is feasible in a 
1.5-dimensional simulation and ensures that the dissipation rate is determined by physical 
effects like a nonlinear cascade and nonlinear steepening rather than by direct numerical 
diffusion (see details in Paper I). 

In order to solve the equations numerically, we use the CIP method (e.g., Yabe, Xiao & 
Utsumi 2001) and the MOCCT method (Stone & Norman 1992). The combination of the 
CIP and MOCCT methods is summarized in Kudoh, Matsumoto & Shibata (1999). The CIP 
method is a useful method to solve an advection equation such as equation (4) accurately, 
and is also applicable to advection terms of hydrodynamic equations. Due to the mirror- 
symmetric boundary condition ai z — 0, the Poisson equation can be simply integrated from 
the midplane of the cloud. 



2.6. Numerical Parameters 

A natural set of fundamental units for this problem are Cso, Hq, and po- These yield a 
time unit to = Hq/csq. The initial magnetic field strength introduces one dimensionless free 
parameter, i.e. 



_ SttPo _ Sttpqc^o 

-Do -Do 



which is the initial ratio of gas to magnetic pressure at z — 0. 

In this cloud, Po is related to the mass-to-flux ratio. For Spitzer's self-gravitating cloud, 
the mass-to-flux ratio normalized to the critical value is 

Hs^2nG'/^^, (23) 
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where 

/oo 
PS dz = 2poHo (24) 
-oo 

is the column density of the Spitzer's self-gravitating cloud. Therefore, 

/3o = fil (25) 

The column density of the cloud we used in this paper is almost equal to that of Spitzer's 
cloud. If we define the column density of the cloud, E, as the integral of density within 
—Zc < z < Zc, then 

E= / p(t = 0)d^~ 0.988 Eg. (26) 

J-Zc 

This means that we can use the value of /is an excellent approximation to the dimensionless 
mass-to-flux ratio of the model cloud. 

The driving force introduces three more free parameters: aa = ad{Ho/c^Q), the dimen- 
sionless amplitude of the acceleration due to driving; i>o = ^^o^o, the dimensionless frequency 
of driving; Za — Za/Ho, the dimensionless scale of the driving region. For simplicity, we take 
Za — 0.1 throughout this paper. 

Dimensional values of all quantities can be found through a choice of Tq and po, along 
with the values of the dimensionless free parameters. For example, if Tq = 10 K and no — 
Po/m = 10^ cm-3, then c^o = 0.2 km s-^ Hq = 0.05 pc, A^5 = Es/m = 3 x lO^^ cm^^ 
to = 2.5 X 10^ yr, and Bq = 20 pG if Po = 1. 



2.7. Oscillatory Motions 

Due to the nonlinearity and inherent randomness in the model, a non-negligible mean 
transverse motion of the cloud can be set up after a few cycles of turbulent driving. In order 
to exclude the energy of the mean transverse cloud motion from our analysis of turbulent 
motions, we introduce (as in Paper I) the oscillating component of the y-velocity, 

v'y^Vy-{Vm(t)), (27) 

where 

Jj^^^^pvydz 

In this equation, Vjn is the mean y— velocity of the cloud and Zf{t) is the full mass position 
of the cloud, which is defined by 

pdz = 0.998^. (29) 



2 



This means that Zf{t) corresponds to the position of the Lagrangian fluid element which is 
initiaUy located at ~ the initial position of the transition region of the temperature. 

In equation (27), we take a time average of Vm for each cycle of the sinusoidal period 
(Tq — I/uo) of the driving force in order to remove the oscillation caused by the driving force. 
The time average, (vm), is deflned in the following manner. For example, (vm) between nTo 
and {n + l)To is calculated as 



where n is an integer. According to the definition, (vm) has the same value between uTq 
and (n + 1)Tq. When the driving force is random, there is really no period. In this case, we 
choose To = to for convenience. 



Table 1 summarizes the parameter study we have performed. The first 5 models with 
different are the same as presented in Paper I. In this paper, we primarily discuss the 
results of varying z/q f3o, as well as the effect of a random driving force. We describe in 
detail a subset of the full range of models listed in Table 1. 



Figure 1 shows the time evolution of models with different frequencies Dq. The left panel 
shows the model with uq = 0.5, the middle panel shows the model with i/q = 1, and the right 
panel shows the model with z/q = 2; all models have 5^ = 30 and /3o = 1- For each model, 
the upper panels show the time evolution of the density. The density profile at various 
times are stacked with time increasing upward in uniform increments of 0.2to- The lower 
panels show the time evolution of v'y at z = 0, where v'y is the oscillating component of the 
y- velocity that we defined in the previous section. As we have shown in Paper I, each density 
evolution plot shows many shock waves propagating in the cloud and significant upward and 
downward motions of the outer portion of the cloud during the period of constant energy 
input {t — lOto — 40fo). Figure 1 shows that the cloud undergoes more expansion and the 
amplitude of v'y becomes greater when the frequency is lower. We return to this point after 
examining the effect of Po- 

Figure 2 presents the same physical quantities as Figure 1, but for models with different 
values of Po- The left panel shows the model with Pq — 0.16, the middle panel shows the 




(30) 



3. 



Results 



3.1. Cloud Dynamics 
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model with /3o = 0.25, and the right panel shows the model with (3o — 4; all models have 
ttfi — 30 and i>o — 1. When /?o = 0.16, the cloud hardly expands, showing no shock waves 
or large up-down motions. The velocity amplitude at 2; = is quite small in comparison 
to other models, and typically subsonic. This model has the strongest magnetic field for a 
given density, and corresponds to the most magnetically subcritical case. When j3o = 0.25, 
the cloud shows large oscillations, although there are not as many shock waves as in the 
standard model with Po — ^ (middle panel of Fig. 1). The case oi Pq — A shows the greatest 
velocity amplitude v'y at the midplane, and many shock waves. This is consistent with the 
trend set by other models that a weaker magnetic field leads to greater disturbances at the 
midplane for a given driving amplitude a^. The (3o — A model is significantly magnetically 
supercritical {ns = 2). However, in contrast to the trend established for increasing Pq by 
the other models, the cloud does not expand nearly as much as in the /?o = 0.25 or /3o = 1 
models. 

The differences in the level of cloud expansion for the models with differing z/q and/or 
(3q can be understood in part by considering the different wavelengths of the input driving. 
The input wavelength Aq is estimated from the Alfven velocity at the midplane Vao and the 
period Tq of the driving force: 



For the models presented in Figure 1 with 9q — 0.5, 1, 2, the corresponding input wavelengths 
are Aq = 2.83ifo, 1.41ifo, O.Tlifo, respectively. The model with uq — 0.5 has an input wave- 
length closest to that of the initial cloud size S.OHq), and this leads to the greatest effect 
of cloud expansion due to the pressure of the nonlinear waves. This correlation is borne out 
further by examining the models of different j3Q presented in Figure 2. When /3o = 0.16 and 
z/q = 1, in addition to the low amplitude fluctuations Vy due to the very strong magnetic 
fleld, we also have Aq ~ 3.5Hq, which is larger than the initial cloud size. Hence, the waves 
are not effective at providing an internal pressure. When (3q = 0.25 and i^o = 1 we find 
Ao ~ 2.83i?o, which is approximately the same as the initial cloud size. On the other hand, 
when Po — A and z>o = 1 we get Aq ~ O.Tlifo, which is significantly smaller than the cloud 
size. Even though there are large velocity fiuctuations Vy at the midplane in this model due 
to the relatively weak magnetic field, the overall expansion of the cloud is inhibited by the 
relatively small ratio Xq/Hq. Our overall results indicate that the cloud suffers less expansion 
when the driving wavelength is larger than the initial cloud size, and the cloud is lifted up 
the most when the driving wavelength is approximately the same as the initial cloud size. 

Figure 3 and Figure 4 show other models which differ from the standard model in the 
values of two separate parameters. The left panel of Figure 3 shows the model with Pq — 0.25 



Ao ~ VaoTq 



Bo 1 V2H^ 







(31) 
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and i>o = 0.5, while the right panel shows the model with /9o = 4 and i>o = 0.5; both models 
have (id = 30. The driving wavelength Aq ~ 5.66i7o for the model in the left panel, and 
Aq ~ IAIHq for the model in the right panel. These two panels are the low frequency 
counterparts of the middle and right panels of Figure 2, respectively. Figure 3 also shows 
that the cloud suffers a lesser expansion when the driving wavelength is larger than the cloud 
size. 

The left panel of Figure 4 shows the model with j3o = 0.25 and Uq = 2, and the right 
panel shows the model with /^o = 4 and Uq = 2; both models have = 30. The driving 
wavelength is Aq ~ 2.83ifo for the model in the left panel, and Aq ~ 0.354iJo foi' the model 
in the right panel. These two models are the high frequency counterparts of the models in 
Figure 3. In contrast to Figure 3, the model with Pq = 0.25 shows greater expansion because 
the wavelength is comparable to the initial cloud size, while the model with (3o — A shows 
very little expansion since the driving wavelength is much smaller than than the initial cloud 
size. 

In Figure 5, we show the results of the models with a random driving force. The 
models have different values of Pq and the same dimensionless amplitude da — 5000. Since 
the random driving contains various frequencies, the time evolution of the density is not 
significantly different from model to model. Each model has a component of the driving at 
the wavelength (which is a function of f3o) that is comparable to the cloud size. However, 
the amplitudes of are a little different in each model. When the magnetic field is strong 
{Po = 0.25), the amplitude is smaller, and vice versa. This result is consistent with the 
results presented in Figure 2. 

To conclude this section, we reemphasize that the wavelength of the driving force is an 
important parameter for the cloud evolution. The driving wavelength which is comparable to 
the cloud size has the greatest effect in expanding the cloud. Thus, the expansion cannot be 
attributed primarily to the commonly used notion of "turbulent pressure", which requires 
that the wavelength of disturbances be much smaller than the typical scale length of the 
cloud (see Bonnazola et al. 1992; Martin, Heyvaerts, & Priest 1998). 



Figure 6 shows the time evolution of energies in the cloud for the case of random driving, 
ad = 5000 and Pq = 1. Each energy is calculated in the following manner: kinetic energy of 
the 2;-component of velocity. 



3.2. Energy Dissipation Rate 




(32) 
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kinetic energy of the y-component of velocity, 



Eky{t) = / -pvldz- Ekm, 



Jo 



(33) 



where Ekm is the kinetic energy of the mean motion of the cloud, i.e.. 



1 /-^/W 
Ekm{t) = -^{vmf / pdz; 



(34) 



^0 



magnetic energy of the y-component of the magnetic field. 



Emit) = [ 



■zj{t) ^2 




(35) 



the sum of the above terms. 



Erit) = Ek,{t) + Eky{t) + Emit). 



(36) 



In the above equations, the integration was done from to Zf{t) because we are interested 
in the energies of the cold material. Strictly speaking, the total energy in each case is twice 
the value we calculate due to the mirror symmetric boundary condition at 2; = 0. Equation 
(33) shows that the mean kinetic energy of the cloud is subtracted from the kinetic energy 
of the the y-component of the velocity, so that E^y is the kinetic energy of the oscillatory 
motions. 

In Figure 6, we examine the evolution of various energies for the random driving force 
model with cid = 5000 and f3o = I- The thick solid line shows Eky, the thin solid line shows 
Em, the dotted line shows E^z, and the dashed-dotted line shows E^- The values of each 
energy are smoothed out over a time t = to to remove small oscillations. Figure 6 is similar to 
the Figure 8 in Paper I, which shows the evolution of these quantities in the standard model 
with sinusoidal driving force. Among the energies, Eky and Em are comparable to each other, 
but Ekz is significantly smaller, i.e., there is approximate equipartition in transverse kinetic 
energy and magnetic fields, and there is much less energy in longitudinal modes. Longitu- 
dinal modes can be generated by the gradient of magnetic field pressure in the ^-direction. 
However, the dominance of the transverse modes implies that the primary dissipation mech- 
anisms arc the nonlinear cascade and nonlinear steepening effects that transfer energy to 
progressively smaller scales until it is dissipated at the grid scale (sec more extensive dis- 
cussion in Paper 1). After the driving force is terminated at t = 40to, the energies decrease 
almost exponentially. The e-folding time is td ~ lOto, which is approximately a sound cross- 
ing time across the cloud size as expanded by the turbulent pressure. The turbulent energy 
is nearly completely dissipated within a few of these crossing times. All of this is consistent 
with the results for sinusoidal driving presented in Paper I. 
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Figure 7 shows the time evolution of Et for different values of i>o with fixed values of 
Po = 1 and ttfi = 30. Figure 8 also shows the time evolution of Et for different values for 
Po but with fixed values z>o = 1 and dd = 30. The energy decreasing time appears to be 
somewhat lesser when the driving frequency is higher. The e-folding time of dissipation for 
each parameter is listed in Table 1. It appears to have a weak correlation with the cloud 
size. The smallest e-folding times {td < 4.0to) occur when the wave length Aq is larger than 
the initial cloud size ~ SHq. In this case, most of the input energy escapes from the cloud 
without affecting it much. Overall, our results show that td ~ lOto is a good approximation 
for most clouds, i.e., that td does not depend very sensitively on most parameters. 



3.3. Spatial Power Spectrum 

Figure 9 shows the spatial power spectra of By and Vy for a model with = 30, Dq = 1, 
and Po = 1. Both are plotted against the dimensionless quantity kHo, where k — 1/i and £ 
is a length scale in the z-direction. The power spectra were taken at t — SOto, in the midst 
of the period of steady turbulent driving, over the spatial range z — — SOHq. There is a 
weak indication of a local peak in the power spectrum of By around A; = I/Aq ~ 0.70/ Hq, 
which corresponds to the driving scale Aq ~ IAIHq. However, the main result here is that 
there is significant power on scales larger than the driving scale {k < 0.70/ Hq) in the spectra 
of both By and Vy. The dotted line shows the dependence k~^^^ for comparison; it is the 
Kolomogorov power spectrum for one- dimensional incompressible hydrodynamic motions 
(e.g., Elmegreen & Scalo 2004). Our simulations of compressible MHD motions shows both 
spectra are similar but have slightly steeper than k"^^^ slopes on small scales(A; > Hq^). 
Figure 10 shows the same power spectra for the random perturbation model with dd — 1000, 
and Po — 1. These power spectra confirm the basic result of the sinusoidal models that there 
is significant power on large scales {k < H^^) and that the slopes are steeper than k~^/^ on 
small scales (A; > Hq^). 

The power spectra show that there is significant power on the largest scales in the cloud 
and that the power spectrum naturally evolves to the form k~^, where p is a positive number, 
even if driving occurs within the cloud on a sub-cloud scale. This large-scale power comes 
from the wave propagation in the gravitationally stratified cloud. When a wave propagates 
into the low density region with greater Alfven speed, the wavelength becomes larger. The 
waves of larger wavelengths also generally have less dissipation. These principles lead to 
significant power on larger scales in a gravitationally stratified cloud. 
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3.4. Correlations between Velocities and Sizes 

Figure 11 plots the time averaged velocity dispersions (o"^)^^ of Lagrangian fluid ele- 
ments within various model clouds as a function of {z)f. where 

^ = ^l\[vl + W+cl (37) 

and {z)t is the time-averaged height of each Lagrangian fluid element. The open circles 
Zp in Table 1) correspond to Lagrangian fluid elements whose initial positions are located at 
z = 2.51Hq, which is close to the edge of the cold cloud. The fllled circles {cth, Zh in Table 
1) correspond to Lagrangian fluid elements whose initial positions are located at z = O.GIHq, 
which is approximately the half-mass position of the cold cloud. An open and filled circle 
is plotted for each model in our complete parameter study; some basic properties of each 
model are summarized in Table 1^. The dotted line shows 

« {z)r- (38) 

This figure shows that the model clouds are in a time-averaged self-gravitational equilibrium 

state under the influence of turbulent pressure (see Paper I). The varying parameters can 
affect the position of the circles, but all model clouds obey the basic correlation shown by the 
dotted line. This is consistent with the observational line-width-size relation of molecular 
clouds (Larson 1981; Myers 1983; Solomon et al. 1987). 

Figure 12 shows the correlation between the velocity dispersion and mean Alfven velocity 
of the cloud Va at the mean position {z)t of the Lagrangian elements. The mean Alfven 
velocity is defined by 

Va 

where 

P 

is the mean density and E is the column density for each Lagrangian element. The dotted 
line shows 

{a')'/' oc Va. (41) 

Results from model clouds with the same are circled and marked. They show that the 
velocity dispersions have a good correlation with the mean Alfven velocity if each cloud has 



^We note that the velocity dispersions listed in Table 1 of Paper I were incorrect, although the Figure 
13 of Paper I was plotted using the correct values of velocity dispersion. This Table contains the correct 
velocity dispersions for all models presented here and in Paper I. 



E 



2(z) 



(39) 
(40) 
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a similar value of /3o, i-e. mass-to-flux ratio. An analysis of the magnetic field measurements 
compiled by Crutcher (1999) shows that there is an excellent correlation between the line- 
of-sight component of the large-scale magnetic field Bi^s and ap^/^ for observed clouds of 
widely varying values of p, Bi^s, and a (Basu 2000). This is essentially the same correlation 
as equation (41); the observational data actually imply an average value a/VA ~ 0.5 (see 
Basu 2005 for discussion). Our models are consistent with this ratio if /3o ~ 1, i-e the mass- 
to-flux ratio is essentially equal to the critical value. A comparison of observed values of 
Bios and column density also reveals that the mass-to-flux ratios are always very close to the 
critical value (Crutcher 1999; Shu et al. 1999; Crutcher 2004; Basu 2005). 



4. Summciry and Discussion 

We have performed 1.5-dimensional numerical simulations of nonlinear MHD waves in 
a gravitationally stratified molecular cloud that is bounded by a hot and tenuous external 
medium. Using the same basic model as presented in Paper I (Kudoh & Basu 2003), we have 
carried out a parameter survey by varying the frequency of the driving force and the magnetic 
field strength of the cloud. We found that the key parameter for the evolution of the cloud 
is the Alfven wavelength of the driving force. If the wavelength is larger than the size of the 
cloud, the cloud is affected less by the waves. The wavelength that is the same order of the 
cloud size is the most effective in expanding the cloud. This means that turbulent expansion 
is different than the usual notion of expansion due to a "turbulent pressure" in which the 
wavelengths need to be much smaller than the cloud size. We also studied the effect of 
a driving force that varies randomly with time and found that significant cloud expansion 
occurs since there is always a component of the driving that is at a favorable frequency for 
expanding the cloud. The evolution of our model clouds reveal the following general features: 
(1) Under the influence of a driving source of Alfvenic disturbances, a cloud shows significant 
upward and downward motions, with an oscillation time scale that is comparable to the cloud 
crossing time; (2) After driving is discontinued, the turbulent energy dissipates exponentially 
with time, with an e-folding time essentially equal to the crossing time across the expanded 
cloud rather than the crossing time of the driving scale; (3) The power spectra show that 
significant power is generated on scales larger than that of the turbulent driving {k < Hq^); 
(4) The hne-width-size relation is obtained by an ensemble of clouds with different physical 
parameters which are individually in a time-averaged self-gravitational equilibrium state. 
The largest amplitude random motions occur in the outer low density regions of a stratified 
cloud. The magnetic field data is best fit by models with /?o ~ 1, which implies that the 
mass-to-flux ratio is very close to the critical value for collapse. 
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In the case of the random driving force, we require a much larger driving amphtude da 
than for the sinusoidal case in order to get a comparable expansion of the cloud (see Table 
1). This is due to the multiplication by a random number in the range [—1, 1] at each time 
step (see eq. [21]) in the random case. The resulting driving efficiency is quite low since 
the driving force F{z,t) is often out of phase with Vy, resulting in a lesser amount of work 
done in the random case than in the sinusoidal case for a given value of a^. A physically 
meaningful comparison of the energetics of the random and sinusoidal cases can be done 
by calculating the energy flux of Alfven waves emerging from just outside the region of the 
artificial driving. We define the Alfven energy fiux in the 2;-direction to be 



and perform the calculation ai z — Za — O.lifo- We find that the time-averaged energy fiux 
of a random driving model with parameters dd = 5000 and /9o = 1 is comparable to (actually 
about twice as large as) that of a sinusoidal driving model with parameters da — 30, (3o — 1, 
and uq — 1. 

For the random driving case (a^ = 5000, — 1), the time-averaged energy fiux mea- 
sured near the surface of the cloud (at the full-mass position) is about 20% of that at z = Za, 
so that about 80% of the input energy is dissipated in the cloud. For the sinusoidal driv- 
ing case {dd = 30, /3o = 1, J^o = 1), the energy flux near the surface is about 30% of that 
at z = Za- The small difference may arise because the random model contains many high 
frequency waves which dissipate quickly. In each case, the dissipated energy flux measured 
while the driving term is applied is comparable to Ex/ta, where Et is the steady-state value 
of the sum of fluctuating magnetic and kinetic energies during the turbulent driving phase, 
and td is the dissipation time measured after turbulent driving is discontinued. Our obtained 
vahic td — lOto, which does not depend sensitively on most parameters, is somewhat longer 
than that estimated from three-dimensional periodic box simulations. We believe that part 
of the discrepancy comes from the generation of longer wavelength modes as the waves travel 
to low-density regions near the cloud's surface. However, one-dimensional simulations are 
also known to have lower dissipation rates than two- or three-dimensional ones (Ostriker et 
al. 2001). Therefore, higher dimensional simulations that include density-stratification of 
the cloud will give the final answer in the future. 

In our ideal MHD model there is no wave damping due to ion-neutral friction, although 
there is wave damping due to a nonlinear cascade and nonlinear steepening, as well as leakage 
of wave energy into the external medium (sec Paper I for a more extensive discussion). We 
expect that ion-neutral friction would be most effective at damping the shortest wavelength 
modes in our simulation (Kulsrud & Pearce 1969), and that the power spectra of By and 
Vy may become steeper at high wavenumber. However, our current results show that the 




(42) 
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wavelengths that are comparable to the initial cloud size are most effective at expanding 
the cloud, so it is not clear exactly how effective ion-neutral friction would be in removing 
turbulent cloud support. Furthermore, we are modeling the large scale evolution of molecular 
clouds with mean column densities of up to a few times 10^^ cm^^, or visual extinction 
Ay ^ 3. At these column densities, the clouds are expected to be significantly ionized by 
background far-ultraviolet starlight so that ion-neutral friction will take a prohibitively long 
time (McKee 1989). Our results show that the dense midplane of the cloud has transonic 
or subsonic motions, and we expect these regions to decouple from the more turbulent 
low density envelope and develop higher column density regions in which the far-ultraviolet 
background is shielded (McKee 1989). In these regions, the ionization fraction is significantly 
lesser, and ambipolar diffusion can enhance the process of core formation. Therefore, models 
of fragmentation of a thin-layer which include the effects of magnetic fields and partial 
ionization due to cosmic rays (Basu & Ciolek 2004) form a complementary approach to ours. 

Our model has the advantage of being global by virtue of including the effect of cloud 
stratification and a cloud boundary, rather than being a local periodic box model. The 
1.5- dimensional approximation has also allowed extremely high resolution which can more 
accurately track the propagation of MHD waves. Our localized turbulent driving source is an 
artificial construct that mimics the various possible sources of disturbances in a real cloud. 
However, a general result is that the wave energy gets quickly redistributed throughout the 
cloud in such a way that most of the energy is on the largest scales. This result provides a 
way out of the "luminosity problem" posed by Basu & Murah (2001), since the dissipation 
rate is then set by the crossing time of the cloud scale rather than that of the internal driving 
scale. Our modeled clouds have a steady-state size ~ 1 pc and large-scale velocity dispersion 
~ 1 km s~^, which makes them comparable to observed dark clouds which can form groupings 
known as the giant molecular clouds (CMC's; see e.g.. Blitz & Williams 1999). The turbulent 
dissipation typically occurs with an e-folding time of a few Myr and is completely dissipated 
by ~ 10^ yr. Even in the case that real molecular clouds and CMC's start their life in a 
turbulent state and then evolve with no further external turbulent input, our results show 
that the turbulent lifetime is less than but within reach of estimated CMC lifetimes, ~ 10^ 
yr; see discussion in Larson (2003). This provides another important motivation for future 
global three-dimensional models which can follow the actual sources of internal turbulent 
driving. As turbulence initially dissipates and a cloud contracts, we expect a transformation 
of some released gravitational potential energy back into MHD wave energy. In addition, 
much of the internal driving may be attributable to the outflows commonly launched from 
the near environs of protostars (see e.g., Bachiller & Tafalla 1999). 
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Fig. 1. — The time evolution for different frequencies (z/q = 0.5, t'o = Ij ^'o = 2) with the 
same parameters of = 30 and f3o = 1. The upper panels show the time evolution of 
density. The density profile at various times are stacked with time increasing upward in 
uniform increments of 0.2io- The lower panels show the time evolution of v'y at z — 0. 




Fig. 2. — The same as Figure 1 but for different plasma beta (/3o = 0.16, (3q = 0.25, /?o = 4) 
with the same dd = 30 and z>o = 1. 
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Fig. 3. — The same as Figure 1 but for different parameters. The left panel shows the case 
of (3q = 0.25, i>o = 0.5 and = 30. The right panel shows the /?o = 4, t>o = 0.5 and dd = 30. 
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Fig. 4. — The same as Figure 1 but for different parameters. The left panel shows the case 
of (5o = 0.25, i>o = 4 and da = 30. The right panel shows the (3o = 4, = 2 and dd = 30. 




Fig. 5. — The time evolution for random driving case with ad = 5000. The left panel shows 
the case of (3^ = 0.25, the middle panel shows /?o = 2, and the right panel shows /3o = 4. 
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Fig. 6. — Time evolution of various energies in the cloud for random driving force with 
= 5000 and jSo = 1. The thick solid line shows E^y, the dotted line shows Ekz, the thin 
solid line shows Em, and dash-dotted line shows Et = Eky + Ekz + Em- The values of each 
energy are smoothed out over a time to in order to remove small oscillations originating from 
the driving force. 
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Fig. 7. — Time evolution of = E^y + E^^ + in the cloud for different frequencies 
(z>o = 0.5, z>o = 1, vq = 2) with the same parameters of hd = 30 and Po = 1. The thick solid 
line shows the case of uq = 0.5, the solid line shows z/q = 1, and the thin solid line shows 
z>o = 2. The values of each energy are smoothed out over every one cycle of the driving force 
in order to remove small oscillations originating from the driving force. 
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Fig. 8. — Time evolution of total energies in the cloud for different plasma beta {(3q = 0.25, 
/?o = 1; /?o = 4) with the same parameters of ad = 30 and z/q = 1. The thick solid line shows 
the case of Pq = 0.25, the solid line shows /3o = 1, and the thin solid line shows Po = 4. The 
values of each energy are smoothed out over every one cycle of the driving force in order to 
remove small oscillations originating from the driving force. 
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Fig. 9. — Spatial power spectra of By and Vy for the case of the sinusoidal driving force with 
parameters = 30, I'o = I and f3o = 1. The spectra were taken from a snapshot at t = 30to- 
The dotted line is proportional to k~^^^ and is shown for comparison. 
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Fig. 10. — Spatial power spectra of By and Vy for the case of the random driving force with 
parameters = 5000 and f3o = 1. The spectra were taken from the snapshot at t — 30to- 
The dotted hne is proportional to k~^^^ and is shown for comparison. 
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Fig. 11. — Time averaged velocity dispersions of different Lagrangian fluid elements for 
different parameters as a function of time averaged positions. The open circles correspond 
to Lagrangian fluid elements whose initial positions are located at z = 2.51Hq, which is close 
to the edge of the cold cloud. The filled circles correspond to Lagrangian fiuid elements whose 
initial positions are located at z = 0.61iJo) which is approximately the half-mass position of 
the cold cloud. The dotted line shows (cr^)j^^ oc {z)f^. Each circle can be associated with a 
particular model in our study by comparison with the numbers in Table 1. 
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Fig. 12. — Time averaged velocity dispersions as a function of the mean Alfven velocity of the 
cloud. The results from the same plasma beta are encircled. The open circles correspond to 
Lagrangian fluid elements whose initial positions are located at z = 2.51Hq. The filled circles 
correspond to Lagrangian fluid elements whose initial positions are located at z = 0.61iJo- 
The dotted line shows {cr'^)]^'^ oc Va- 



